EMA Modeling

load packages

library(tidyverse)
library(lmerTest)

define palettes

pal_process = wesanderson::wes_palette("Zissou1", 6, "continuous")

source data

source("load_data.R")

standardize brain data

# standardize and winsorize
betas_std = betas %>%
  group_by(roi, session) %>%
  mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
         meanPE_std = ifelse(meanPE_std > 3, 3,
                      ifelse(meanPE_std < -3, -3, meanPE_std))) %>%
  ungroup()

dots_std = dots %>%
  group_by(map, test, mask, session) %>%
  mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
         dotProduct_std = ifelse(dotProduct_std > 3, 3,
                          ifelse(dotProduct_std < -3, -3, dotProduct_std))) %>%
  ungroup()

# join data frames
dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
  mutate(subjectID = as.character(subjectID)) %>%
  ungroup() %>%
  mutate(subjectID = as.factor(subjectID),
         condition = as.factor(condition),
         control = as.factor(control),
         roi = as.factor(roi),
         process = as.factor(process),
         test = as.factor(test))

tidy data

ema_tidy = ema %>%
  select(subjectID, desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m, hunger.m) %>%
  unique()

model_data = dataset %>%
  filter(((!grepl("neuralsig", map) & test == "association" & mask == "masked") | (grepl("neuralsig", map) & mask == "unmasked")) & 
           session == "all" & control == "nature" & condition == "snack") %>%
  unique() %>%
  group_by(process, subjectID) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
         processPE = paste0(process, "_ROI"),
         processPEV = paste0(process, "_PEV")) %>%
  ungroup() %>%
  select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map, process)) %>%
  spread(processPEV, dotProduct_std) %>%
  spread(processPE, meanProcessPEstd) %>%
  group_by(subjectID) %>%
  fill(everything(), .direction = "down") %>%
  fill(everything(), .direction = "up") %>%
  select(-c(roi, craving_ROI, craving_regulation_ROI, test, mask)) %>%
  unique() %>%
  mutate(balance_ROI = cognitive_control_ROI - reward_ROI) %>%
  right_join(., ema_tidy, by = "subjectID") %>%
  na.omit()

visualize

ROIs

model_data %>%
  select(subjectID, contains("ROI"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
  unique() %>%
  gather(neuro_var, neuro_val, contains("ROI")) %>%
  mutate(neuro_var = factor(neuro_var, levels = c("balance_ROI", "cognitive_control_ROI",
                                                  "reward_ROI", "value_ROI"))) %>%
  gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
  group_by(ema_var) %>%
  do({
    plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
      geom_point(alpha = .5) +
      geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
      scale_color_manual(name = "", values = pal_process) +
      labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
      theme_minimal(base_size = 14)
    
    print(plot)
    data.frame()
  })

MAMs

model_data %>%
  select(subjectID, contains("PEV"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
  unique() %>%
  gather(neuro_var, neuro_val, contains("PEV")) %>%
  mutate(neuro_var = factor(neuro_var, levels = c("cognitive_control_PEV", "reward_PEV", "value_PEV",
                                                  "craving_PEV", "craving_regulation_PEV"))) %>%
  gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
  group_by(ema_var) %>%
  do({
    plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
      geom_point(alpha = .5) +
      geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
      scale_color_manual(name = "", values = pal_process[2:6]) +
      labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
      theme_minimal(base_size = 14)
    
    print(plot)
    data.frame()
  })

run regression models

specify model variables

options(na.action = "na.fail")
data_ctrl = caret::trainControl(method = "repeatedcv", number = 5, repeats = 5)

desire presence proportion

ROI models

desire_prop_roi_full = lm(desire_pres0.prop ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data)
(desire_prop_roi_mods = MuMIn::dredge(desire_prop_roi_full))

MuMIn::get.models(desire_prop_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_desire_prop_roi = caret::train(desire_pres0.prop ~ reward_ROI,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

MAM models

desire_prop_mam_full = lm(desire_pres0.prop ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data)
(desire_prop_mam_mods = MuMIn::dredge(desire_prop_mam_full))

MuMIn::get.models(desire_prop_mam_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_desire_prop_mam = caret::train(desire_pres0.prop ~ value_PEV,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

desire_prop_null = lm(desire_pres0.prop ~ 1, data = model_data)

desire_prop_combined = caret::train(desire_pres0.prop ~ reward_ROI + value_PEV,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

AIC(desire_prop_null, best_desire_prop_roi$finalModel, best_desire_prop_mam$finalModel, desire_prop_combined$finalModel) %>%
  rownames_to_column() %>%
  arrange(AIC)

Summarize the overall best fitting

summary(best_desire_prop_mam$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27760 -0.12561 -0.01319  0.10406  0.37274 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33844    0.01682  20.119   <2e-16 ***
## value_PEV    0.02656    0.01382   1.922   0.0576 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1494 on 95 degrees of freedom
## Multiple R-squared:  0.03743,    Adjusted R-squared:  0.0273 
## F-statistic: 3.694 on 1 and 95 DF,  p-value: 0.05761

best_desire_prop_mam

## Linear Regression 
## 
## 97 samples
##  1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 77, 78, 78, 78, 77, 78, ... 
## Resampling results:
## 
##   RMSE       Rsquared    MAE      
##   0.1496044  0.09099185  0.1255609
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(best_desire_prop_mam$resample$Rsquared)

## [1] 0.102751

sd(best_desire_prop_mam$resample$RMSE)

## [1] 0.01693802

enactment proportion

ROI models

enact_prop_roi_full = lm(enact_prop ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data)
(enact_prop_roi_mods = MuMIn::dredge(enact_prop_roi_full))

MuMIn::get.models(enact_prop_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_enact_prop_roi = caret::train(enact_prop ~ reward_ROI,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

MAM models

enact_prop_mam_full = lm(enact_prop ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data)
(enact_prop_mam_mods = MuMIn::dredge(enact_prop_mam_full))

MuMIn::get.models(enact_prop_mam_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_enact_prop_mam = caret::train(enact_prop ~ cognitive_control_PEV + craving_regulation_PEV + reward_PEV + value_PEV,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

enact_prop_null = lm(enact_prop ~ 1, data = model_data)

enact_prop_combined = caret::train(enact_prop ~ reward_ROI + cognitive_control_PEV + craving_regulation_PEV + reward_PEV + value_PEV,
                     data = model_data,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

AIC(enact_prop_null, best_enact_prop_roi$finalModel, best_enact_prop_mam$finalModel, enact_prop_combined$finalModel) %>%
  rownames_to_column() %>%
  arrange(AIC)

Summarize the overall best fitting

summary(best_enact_prop_mam$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.61127 -0.13058  0.00977  0.10827  0.49466 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.47037    0.02893  16.261  < 2e-16 ***
## cognitive_control_PEV  -0.06868    0.02749  -2.498  0.01426 *  
## craving_regulation_PEV  0.03779    0.02324   1.626  0.10731    
## reward_PEV              0.19511    0.06939   2.812  0.00602 ** 
## value_PEV              -0.12959    0.06221  -2.083  0.04001 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2121 on 92 degrees of freedom
## Multiple R-squared:  0.1349, Adjusted R-squared:  0.09732 
## F-statistic: 3.587 on 4 and 92 DF,  p-value: 0.009164

best_enact_prop_mam

## Linear Regression 
## 
## 97 samples
##  4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 78, 77, 78, 77, 78, 78, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE      
##   0.214186  0.1087687  0.1628077
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(best_enact_prop_mam$resample$Rsquared)

## [1] 0.102848

sd(best_enact_prop_mam$resample$RMSE)

## [1] 0.03547659